Ejemplo 6
Inicialización
Datos globales para modelos
#--- Ejemplo 6 ---
#-Reading data-
desastres<-read.table("desastres.txt",header=TRUE)
n<-nrow(desastres)
plot(desastres,type="l")
plot(desastres[2:n,2]-desastres[1:(n-1),2],type="l")
plot(log(desastres[2:n,2])-log(desastres[1:(n-1),2]),type="l")
#-Defining data-
data<-list("n"=n,"y"=desastres$No.Desastres,"x"=desastres$Anho)
data<-list("n"=n,"y"=c(desastres$No.Desastres[1:(n-6)],rep(NA,6)),"x"=desastres$Anho)
Funciones de graficación
source("util.R")
graficarY1 <- function(modelo){
if("BUGSoutput" %in% names(modelo)){
modelo <- modelo$BUGSoutput
}
out.sum<-modelo$summary
#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(desastres[,2],out.yf[,c(1,3,7)])
ymax<-max(desastres[,2],out.yf[,c(1,3,7)])
par(mfrow=c(1,1))
plot(desastres,type="l",col="grey80",ylim=c(ymin,ymax))
lines(desastres[,1],out.yf[,1],lwd=2,col=2)
lines(desastres[,1],out.yf[,3],lty=2,col=2)
lines(desastres[,1],out.yf[,7],lty=2,col=2)
lines(desastres[,1],out.yf[,5],lwd=2,col=4)
#Medias
out.mu<-out.sum[grep("mu",rownames(out.sum)),]
par(mfrow=c(1,1))
plot(desastres,type="l",col="grey80")
lines(desastres[,1],out.mu[,1],lwd=2,col=2)
lines(desastres[,1],out.mu[,3],lty=2,col=2)
lines(desastres[,1],out.mu[,7],lty=2,col=2)
}
graficarTasa <- function(modelo){
if("BUGSoutput" %in% names(modelo)){
modelo <- modelo$BUGSoutput
}
out.sum<-modelo$summary
out.tasa<-out.sum[grep("lambda",rownames(out.sum)),]
out.tasa<-out.sum[grep("p",rownames(out.sum)),]
or<-order(mortality$x)
ymin<-min(mortality$y/mortality$n,out.tasa[,c(1,3,7)])
ymax<-max(mortality$y/mortality$n,out.tasa[,c(1,3,7)])
plot(mortality$x,mortality$y/mortality$n,ylim=c(ymin,ymax))
lines(mortality$x[or],out.tasa[or,1],lwd=2,col=4)
lines(mortality$x[or],out.tasa[or,3],lty=2,col=4)
lines(mortality$x[or],out.tasa[or,7],lty=2,col=4)
}
Ejercicio 6a
cat Ej6a.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## # y[i] ~ dnegbin(p[i],r)
## # logit(p[i])<-beta[1]+beta[2]*x[i]
## # mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## y[i] ~ dpois(mu[i])
## log(mu[i])<-beta[1]+beta[2]*x[i]
## }
## #Priors
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## #aux ~ dpois(5)
## #r <- aux + 1
##
## #Prediction 1
## #Neg Binomial
## #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6a.txt",
n.iter=20000,n.chains=1,n.burnin=4000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6a.txt",
n.iter=20000,n.chains=4,n.burnin=4000,n.thin=5)
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 106
## Unobserved stochastic nodes: 120
## Total graph size: 679
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)
Betas
betasMCMC(ej6a.bugs)
## mean 2.5% 97.5% prob
## beta[1] 34.97886250 25.40000 44.99000 0
## beta[2] -0.01816264 -0.02347 -0.01311 0
betasMCMC(ej6a.jags)
## mean 2.5% 97.5% prob
## beta[1] 35.04765163 25.49703245 44.9492355 0
## beta[2] -0.01819738 -0.02343192 -0.0131291 0
Ejercicio 6a - Predicción
cat Ej6a.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## # y[i] ~ dnegbin(p[i],r)
## # logit(p[i])<-beta[1]+beta[2]*x[i]
## # mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## y[i] ~ dpois(mu[i])
## log(mu[i])<-beta[1]+beta[2]*x[i]
## }
## #Priors
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## #aux ~ dpois(5)
## #r <- aux + 1
##
## #Prediction 1
## #Neg Binomial
## #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
data<-list("n"=n,"y"=c(desastres$No.Desastres[1:(n-6)],rep(NA,6)),"x"=desastres$Anho)
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6a.txt",
n.iter=20000,n.chains=1,n.burnin=4000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6a.txt",
n.iter=20000,n.chains=4,n.burnin=4000,n.thin=5)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 106
## Unobserved stochastic nodes: 120
## Total graph size: 679
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)
Betas
betasMCMC(ej6a.bugs)
## mean 2.5% 97.5% prob
## beta[1] 34.97886250 25.40000 44.99000 0
## beta[2] -0.01816264 -0.02347 -0.01311 0
betasMCMC(ej6a.jags)
## mean 2.5% 97.5% prob
## beta[1] 35.05034708 25.64597116 45.11129741 7.8125e-05
## beta[2] -0.01819942 -0.02351621 -0.01321869 7.8125e-05
DIC
dicMCMC(ej6a.bugs)
## [1] 338.1
dicMCMC(ej6a.jags)
## [1] 339.3104
Restore data
#reestablecer para no tener predicción
data<-list("n"=n,"y"=desastres$No.Desastres,"x"=desastres$Anho)
Ejercicio 6aa
cat Ej6aa.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## y[i] ~ dnegbin(p[i],r)
## logit(p[i])<-beta[1]+beta[2]*x[i]
## mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## # y[i] ~ dpois(mu[i])
## # log(mu[i])<-beta[1]+beta[2]*x[i]
## }
## #Priors
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## aux ~ dpois(10)
## r <- aux + 1
##
## #Prediction 1
## #Neg Binomial
## for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## #for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6aa.txt",
n.iter=10000,n.chains=2,n.burnin=1000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6aa.txt",
n.iter=10000,n.chains=2,n.burnin=1000,n.thin=1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 112
## Unobserved stochastic nodes: 115
## Total graph size: 1244
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
Gráficas
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)
exploracionMCMC_r(ej6a.bugs)
exploracionMCMC_r(ej6a.jags)
Betas
betasMCMC(ej6a.bugs)
## mean 2.5% 97.5% prob
## beta[1] -33.72750056 -43.80000 -23.9200 0
## beta[2] 0.01878268 0.01356 0.0241 0
betasMCMC(ej6a.jags)
## mean 2.5% 97.5% prob
## beta[1] -18.04765706 -25.63400634 -5.63592620 0
## beta[2] 0.01047011 0.00391693 0.01451908 0
Ejercicio 6b
cat Ej6b.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## # y[i] ~ dnegbin(p[i],r)
## # logit(p[i])<-beta[1]+beta[2]*step(x[i]-tau)
## # mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## y[i] ~ dpois(mu[i])
## log(mu[i])<-beta[1]+beta[2]*step(x[i]-tau)
## }
## #Priors
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## #aux ~ dpois(5)
## #r <- aux + 1
## aux2 ~ dcat(a[])
## tau <- aux2 + 1850
## for (j in 1:112) { a[j]<- 1/112}
## #Prediction 1
## #Neg Binomial
## #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6b.bugs<-bugs(data,inits,parameters,model.file="Ej6b.txt",
n.iter=5000,n.chains=1,n.burnin=500)
ej6b.jags<-jags(data,inits,parameters,model.file="Ej6b.txt",
n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 112
## Unobserved stochastic nodes: 115
## Total graph size: 1133
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6b.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6b.jags)
Betas
betasMCMC(ej6b.bugs)
## mean 2.5% 97.5% prob
## beta[1] 1.142430 0.951585 1.320000 0
## beta[2] -1.261346 -1.570525 -0.965045 0
betasMCMC(ej6b.jags)
## mean 2.5% 97.5% prob
## beta[1] 1.078548 0.9452293 1.3186034 0.005555556
## beta[2] -1.170457 -1.5511344 -0.9268255 0.005555556
Ejercicio 6bb
cat Ej6bb.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## y[i] ~ dnegbin(p[i],r)
## logit(p[i])<-beta[1]+beta[2]*step(x[i]-tau)
## mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## # y[i] ~ dpois(mu[i])
## # log(mu[i])<-beta[1]+beta[2]*step(x[i]-tau)
## }
## #Priors
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## aux ~ dpois(5)
## r <- aux + 1
## aux2 ~ dcat(a[])
## tau <- aux2 + 1850
## for (j in 1:112) { a[j]<- 1/112}
## #Prediction 1
## #Neg Binomial
## for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## #for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6b.bugs<-bugs(data,inits,parameters,model.file="Ej6bb.txt",
n.iter=5000,n.chains=2,n.burnin=1000)
ej6b.jags<-jags(data,inits,parameters,model.file="Ej6bb.txt",
n.iter=5000,n.chains=2,n.burnin=1000,n.thin=5)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 112
## Unobserved stochastic nodes: 116
## Total graph size: 1810
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6b.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6b.jags)
Betas
betasMCMC(ej6b.bugs)
## mean 2.5% 97.5% prob
## beta[1] 0.8680623 0.2378725 1.415025 0.00425
## beta[2] 1.2579773 0.9203925 1.590000 0.00000
betasMCMC(ej6b.jags)
## mean 2.5% 97.5% prob
## beta[1] 0.8469051 0.2165423 1.427110 0.00375
## beta[2] 1.2590502 0.9339628 1.603908 0.00000
Ejercicio 6c
cat Ej6c.txt
## model
## {
## #Likelihood
## #Space eq.
## for (i in 1:n) {
## y[i] ~ dpois(mu[i])
## log(mu[i])<-beta[i]
## }
## #State eq.
## for (i in 2:n) {
## beta[i] ~ dnorm(beta[i-1],tau.b)
## }
## #Priors
## beta[1] ~ dnorm(0,0.001)
## tau.b ~ dgamma(0.001,0.001)
##
## #Prediction 1
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
# inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
parameters<-c("beta", "tau.b","yf1","mu")
ej6c.bugs<-bugs(data,inits,parameters,model.file="Ej6c.txt",
n.iter=5000,n.chains=1,n.burnin=500)
ej6c.jags<-jags(data,inits,parameters,model.file="Ej6c.txt",
n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "x" in data
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 112
## Unobserved stochastic nodes: 225
## Total graph size: 454
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6c.bugs)
print("Jags")
## [1] "Jags"
exploracionMCMC(ej6c.jags)
Betas
betasMCMC2 <- function(modelo){
if("BUGSoutput" %in% names(modelo)){
modelo <- modelo$BUGSoutput
}
out <- modelo$sims.list
out.sum<-modelo$summary
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
dimnames(out.sum.t)[[2]][4]<-"prob"
# print(modelo$sims.list$beta)
print(out.sum.t)
}
betasMCMC2(ej6c.bugs)
## mean 2.5% 97.5% prob
## beta[1] 1.186404476 0.64511750 1.71700000 0.0004444444
## beta[2] 1.176569556 0.67658493 1.66252492 0.0000000000
## beta[3] 1.124161800 0.64859499 1.57500000 0.0000000000
## beta[4] 1.049636000 0.56629499 1.48052492 0.0000000000
## beta[5] 1.019035578 0.54708000 1.44000000 0.0002222222
## beta[6] 1.055533889 0.58808997 1.48352492 0.0000000000
## beta[7] 1.066987567 0.62878493 1.48852492 0.0000000000
## beta[8] 1.077790578 0.63284248 1.49000000 0.0000000000
## beta[9] 1.064095044 0.60964750 1.47505000 0.0002222222
## beta[10] 1.120860822 0.68908997 1.54100000 0.0000000000
## beta[11] 1.114888267 0.68958493 1.53700000 0.0000000000
## beta[12] 1.111782956 0.69667482 1.49452492 0.0000000000
## beta[13] 1.112582609 0.68458997 1.50500000 0.0000000000
## beta[14] 1.087464256 0.66088250 1.48052500 0.0004444444
## beta[15] 1.134025000 0.71293746 1.52152492 0.0000000000
## beta[16] 1.206312933 0.79754249 1.61300000 0.0000000000
## beta[17] 1.213023956 0.80384249 1.62500000 0.0000000000
## beta[18] 1.230625867 0.83744428 1.62704969 0.0000000000
## beta[19] 1.261067356 0.86422233 1.66252492 0.0000000000
## beta[20] 1.256854356 0.87554750 1.67852493 0.0000000000
## beta[21] 1.242330578 0.85782738 1.65852492 0.0000000000
## beta[22] 1.197264644 0.79614249 1.59000000 0.0000000000
## beta[23] 1.162509911 0.75122230 1.54752492 0.0000000000
## beta[24] 1.178227378 0.76949499 1.57900000 0.0000000000
## beta[25] 1.179947733 0.77742231 1.56900000 0.0000000000
## beta[26] 1.164413933 0.75654750 1.55100000 0.0000000000
## beta[27] 1.203783556 0.79613746 1.61852492 0.0000000000
## beta[28] 1.206855889 0.79754937 1.62652492 0.0000000000
## beta[29] 1.172528311 0.76814750 1.59652492 0.0000000000
## beta[30] 1.145873956 0.73849499 1.58000000 0.0000000000
## beta[31] 1.097180844 0.67338493 1.49852492 0.0000000000
## beta[32] 1.078060978 0.66419499 1.48757425 0.0000000000
## beta[33] 1.004883378 0.57760000 1.41452491 0.0000000000
## beta[34] 0.950060960 0.51099499 1.36552491 0.0000000000
## beta[35] 0.911446511 0.46928500 1.34352500 0.0002222222
## beta[36] 0.865600564 0.41949500 1.30705000 0.0002222222
## beta[37] 0.781966251 0.33398000 1.22000000 0.0013333333
## beta[38] 0.700140861 0.23969000 1.11905000 0.0015555556
## beta[39] 0.641424687 0.18442250 1.08452500 0.0037777778
## beta[40] 0.557892641 0.06231775 0.99670500 0.0157777778
## beta[41] 0.467437037 -0.05962300 0.92395250 0.0377777778
## beta[42] 0.393509395 -0.17106750 0.86500750 0.0711111111
## beta[43] 0.332014646 -0.22832000 0.80730500 0.1057777778
## beta[44] 0.279232468 -0.27471500 0.77456000 0.1437777778
## beta[45] 0.233980967 -0.32697750 0.73136750 0.1811111111
## beta[46] 0.197431397 -0.35985250 0.69696250 0.2191111111
## beta[47] 0.118350283 -0.44276250 0.61862500 0.3093333333
## beta[48] 0.063620082 -0.52201250 0.57589750 0.3786666667
## beta[49] 0.037567481 -0.54387750 0.55365250 0.4197777778
## beta[50] 0.007925561 -0.57476250 0.52636750 0.4682222222
## beta[51] -0.005325450 -0.61120500 0.52350000 0.4820000000
## beta[52] -0.017300148 -0.63160000 0.53476750 0.4902222222
## beta[53] -0.021949274 -0.63070500 0.53286250 0.4733333333
## beta[54] -0.002374363 -0.57310500 0.55223000 0.4964444444
## beta[55] 0.039074907 -0.55871000 0.61164500 0.4484444444
## beta[56] 0.031045862 -0.54865000 0.59510000 0.4577777778
## beta[57] 0.032573767 -0.53162000 0.59450000 0.4413333333
## beta[58] 0.060416980 -0.48085250 0.61166000 0.4064444444
## beta[59] 0.041604198 -0.51738250 0.57896750 0.4340000000
## beta[60] 0.003886435 -0.55991000 0.54470000 0.4828888889
## beta[61] -0.055997085 -0.61332500 0.46666250 0.4322222222
## beta[62] -0.093295581 -0.64683500 0.42205750 0.3753333333
## beta[63] -0.134151626 -0.69679250 0.37190500 0.3017777778
## beta[64] -0.176386264 -0.75173500 0.34095250 0.2568888889
## beta[65] -0.222535966 -0.79571000 0.29851000 0.2191111111
## beta[66] -0.250927162 -0.85855250 0.27646750 0.1902222222
## beta[67] -0.283007907 -0.89337250 0.26708750 0.1617777778
## beta[68] -0.295794442 -0.89886750 0.25710500 0.1375555556
## beta[69] -0.313690998 -0.91796250 0.19672000 0.1228888889
## beta[70] -0.310900799 -0.91280500 0.20661000 0.1246666667
## beta[71] -0.284486115 -0.87571000 0.21991000 0.1457777778
## beta[72] -0.243625697 -0.82836250 0.24407250 0.1822222222
## beta[73] -0.235973155 -0.85269750 0.25327750 0.1908888889
## beta[74] -0.233601357 -0.85687500 0.26327250 0.1900000000
## beta[75] -0.216706930 -0.87686750 0.26910500 0.2146666667
## beta[76] -0.180229579 -0.79441500 0.29446250 0.2617777778
## beta[77] -0.125199487 -0.71162750 0.37233250 0.3275555556
## beta[78] -0.073217687 -0.63300500 0.45077750 0.3915555556
## beta[79] -0.025089270 -0.56238750 0.51972000 0.4540000000
## beta[80] 0.055059698 -0.46825250 0.61167250 0.4246666667
## beta[81] 0.113543070 -0.39249750 0.66156250 0.3391111111
## beta[82] 0.153783501 -0.34525250 0.70487750 0.2891111111
## beta[83] 0.147676497 -0.34041000 0.69402500 0.3040000000
## beta[84] 0.150433047 -0.34150250 0.69286250 0.3068888889
## beta[85] 0.155751608 -0.34341500 0.70704500 0.2977777778
## beta[86] 0.146421960 -0.36305750 0.70730500 0.3237777778
## beta[87] 0.141963327 -0.35497250 0.69420500 0.3168888889
## beta[88] 0.145182792 -0.33694000 0.68863500 0.3162222222
## beta[89] 0.154931453 -0.33313500 0.70361000 0.3066666667
## beta[90] 0.168015509 -0.35680500 0.73011000 0.2915555556
## beta[91] 0.162448044 -0.35905750 0.74795000 0.3004444444
## beta[92] 0.088263030 -0.44461500 0.64138750 0.3833333333
## beta[93] -0.012043187 -0.56198250 0.52694500 0.4815555556
## beta[94] -0.089814005 -0.66776250 0.45510500 0.3864444444
## beta[95] -0.147530418 -0.74976250 0.39850000 0.3208888889
## beta[96] -0.181250019 -0.82110000 0.37100750 0.2728888889
## beta[97] -0.213382314 -0.85020000 0.33967250 0.2315555556
## beta[98] -0.324215623 -1.02105000 0.23230000 0.1326666667
## beta[99] -0.419107433 -1.15657500 0.16766250 0.0835555556
## beta[100] -0.499654283 -1.30157500 0.10521000 0.0533333333
## beta[101] -0.558980694 -1.40857500 0.04982850 0.0404444444
## beta[102] -0.630100846 -1.47552500 -0.00975905 0.0226666667
## beta[103] -0.685981070 -1.55352500 -0.03760975 0.0184444444
## beta[104] -0.733273178 -1.64052500 -0.05942825 0.0137777778
## beta[105] -0.765173364 -1.67500000 -0.11404250 0.0084444444
## beta[106] -0.785009173 -1.76757500 -0.13663750 0.0062222222
## beta[107] -0.795944328 -1.80200000 -0.13360000 0.0051111111
## beta[108] -0.822018586 -1.87605000 -0.15024750 0.0060000000
## beta[109] -0.845735502 -1.92062500 -0.13918750 0.0051111111
## beta[110] -0.856390510 -1.98262500 -0.10969750 0.0077777778
## beta[111] -0.878551151 -2.08205000 -0.10609000 0.0091111111
## beta[112] -0.888979448 -2.11862500 -0.08120125 0.0142222222
betasMCMC2(ej6c.jags)
## mean 2.5% 97.5% prob
## beta[1] 1.163482e+00 0.61653332 1.688036414 0.000000000
## beta[2] 1.155103e+00 0.69062415 1.661144926 0.000000000
## beta[3] 1.104094e+00 0.62047084 1.542039619 0.000000000
## beta[4] 1.036272e+00 0.55949126 1.486278335 0.000000000
## beta[5] 1.012751e+00 0.49804356 1.459594189 0.001111111
## beta[6] 1.059754e+00 0.62721248 1.472797977 0.000000000
## beta[7] 1.084537e+00 0.69323320 1.508018336 0.000000000
## beta[8] 1.098809e+00 0.69791268 1.486678246 0.000000000
## beta[9] 1.086438e+00 0.67627994 1.512372381 0.000000000
## beta[10] 1.144797e+00 0.75348161 1.535781203 0.000000000
## beta[11] 1.129157e+00 0.73143817 1.525132094 0.000000000
## beta[12] 1.120300e+00 0.70079123 1.504650873 0.000000000
## beta[13] 1.114881e+00 0.70143074 1.497270777 0.000000000
## beta[14] 1.073313e+00 0.64792592 1.438048293 0.000000000
## beta[15] 1.117922e+00 0.69650869 1.481535467 0.000000000
## beta[16] 1.189446e+00 0.79861095 1.590036396 0.000000000
## beta[17] 1.185427e+00 0.73011400 1.598368200 0.000000000
## beta[18] 1.203445e+00 0.81508297 1.585723633 0.000000000
## beta[19] 1.240088e+00 0.81907015 1.697939469 0.000000000
## beta[20] 1.235901e+00 0.82777764 1.679552441 0.000000000
## beta[21] 1.217214e+00 0.79915733 1.616699970 0.000000000
## beta[22] 1.171920e+00 0.78091496 1.533623247 0.000000000
## beta[23] 1.136336e+00 0.71349773 1.536792642 0.000000000
## beta[24] 1.161959e+00 0.75188490 1.575042077 0.000000000
## beta[25] 1.167846e+00 0.76108877 1.574735262 0.000000000
## beta[26] 1.150074e+00 0.75437353 1.545563593 0.000000000
## beta[27] 1.192170e+00 0.81183412 1.596889877 0.000000000
## beta[28] 1.191645e+00 0.83814811 1.588207551 0.000000000
## beta[29] 1.164381e+00 0.76911834 1.592200963 0.000000000
## beta[30] 1.129464e+00 0.77283331 1.528607882 0.000000000
## beta[31] 1.080048e+00 0.70840429 1.464420587 0.000000000
## beta[32] 1.051264e+00 0.66552333 1.457800381 0.000000000
## beta[33] 9.813956e-01 0.59429529 1.396188831 0.000000000
## beta[34] 9.266961e-01 0.51515735 1.329929522 0.000000000
## beta[35] 8.829496e-01 0.44839284 1.303133477 0.000000000
## beta[36] 8.359225e-01 0.41095324 1.294551129 0.000000000
## beta[37] 7.492307e-01 0.33146463 1.157445192 0.000000000
## beta[38] 6.740897e-01 0.23423120 1.110356486 0.002222222
## beta[39] 6.145917e-01 0.17257305 1.021404194 0.003333333
## beta[40] 5.361867e-01 0.05319936 0.943473359 0.013333333
## beta[41] 4.498741e-01 -0.03710060 0.875984208 0.042222222
## beta[42] 3.819056e-01 -0.15427053 0.802354605 0.073333333
## beta[43] 3.259579e-01 -0.18514443 0.795971909 0.088888889
## beta[44] 2.815060e-01 -0.21435915 0.724862365 0.110000000
## beta[45] 2.330995e-01 -0.26295307 0.657488049 0.152222222
## beta[46] 1.921749e-01 -0.30005829 0.630028311 0.217777778
## beta[47] 1.136508e-01 -0.39280931 0.566258514 0.307777778
## beta[48] 5.962312e-02 -0.48034487 0.533535380 0.394444444
## beta[49] 3.093576e-02 -0.47426470 0.501936724 0.423333333
## beta[50] -6.078716e-04 -0.57446119 0.485795194 0.478888889
## beta[51] -4.092496e-03 -0.60872509 0.483078815 0.491111111
## beta[52] -2.026932e-03 -0.56205485 0.508664677 0.477777778
## beta[53] 4.010552e-05 -0.59312281 0.521204908 0.475555556
## beta[54] 2.270509e-02 -0.56487598 0.528294140 0.465555556
## beta[55] 7.644045e-02 -0.45451180 0.572506719 0.365555556
## beta[56] 7.692947e-02 -0.43879035 0.567063112 0.368888889
## beta[57] 8.484146e-02 -0.44327232 0.565981430 0.348888889
## beta[58] 1.180576e-01 -0.41768240 0.578534875 0.294444444
## beta[59] 9.917300e-02 -0.44412396 0.589258202 0.310000000
## beta[60] 4.869441e-02 -0.50816172 0.563041565 0.386666667
## beta[61] -1.430205e-02 -0.62021154 0.497161969 0.472222222
## beta[62] -5.177988e-02 -0.65583266 0.438087809 0.461111111
## beta[63] -1.038017e-01 -0.71483654 0.414743254 0.367777778
## beta[64] -1.500083e-01 -0.74871144 0.351465317 0.301111111
## beta[65] -2.065042e-01 -0.78662445 0.338876210 0.236666667
## beta[66] -2.299535e-01 -0.86039571 0.298419469 0.222222222
## beta[67] -2.790901e-01 -0.93307976 0.262859338 0.196666667
## beta[68] -2.938798e-01 -0.96122647 0.254392398 0.184444444
## beta[69] -3.174843e-01 -0.97663053 0.246526410 0.150000000
## beta[70] -3.231845e-01 -1.01894244 0.248645459 0.152222222
## beta[71] -3.074042e-01 -0.93658145 0.231949022 0.151111111
## beta[72] -2.691884e-01 -0.90162595 0.267300871 0.190000000
## beta[73] -2.629401e-01 -0.90989779 0.260624238 0.185555556
## beta[74] -2.656030e-01 -0.92635537 0.232333012 0.190000000
## beta[75] -2.535037e-01 -0.90792325 0.249824727 0.197777778
## beta[76] -2.155179e-01 -0.82136378 0.291537114 0.244444444
## beta[77] -1.542421e-01 -0.76890759 0.323221445 0.313333333
## beta[78] -9.736057e-02 -0.68082126 0.401125747 0.377777778
## beta[79] -4.160857e-02 -0.57936692 0.467998791 0.452222222
## beta[80] 3.918525e-02 -0.47499596 0.585720105 0.427777778
## beta[81] 1.018864e-01 -0.45137656 0.639037576 0.355555556
## beta[82] 1.453024e-01 -0.39149125 0.714086048 0.297777778
## beta[83] 1.382575e-01 -0.42858932 0.694016493 0.308888889
## beta[84] 1.306509e-01 -0.36596773 0.686906374 0.333333333
## beta[85] 1.433223e-01 -0.35410625 0.684472438 0.307777778
## beta[86] 1.383640e-01 -0.35133408 0.683252307 0.314444444
## beta[87] 1.375278e-01 -0.36178998 0.698959365 0.310000000
## beta[88] 1.437146e-01 -0.38701425 0.729229515 0.298888889
## beta[89] 1.511331e-01 -0.37243699 0.713111338 0.291111111
## beta[90] 1.754862e-01 -0.33596744 0.758854117 0.270000000
## beta[91] 1.732897e-01 -0.34962951 0.730623203 0.288888889
## beta[92] 1.142001e-01 -0.43139141 0.676031217 0.355555556
## beta[93] 1.366668e-02 -0.49220838 0.567634432 0.476666667
## beta[94] -6.193843e-02 -0.58572179 0.464664670 0.417777778
## beta[95] -1.060191e-01 -0.62126391 0.387466557 0.354444444
## beta[96] -1.321277e-01 -0.67082218 0.360501039 0.340000000
## beta[97] -1.555074e-01 -0.71298970 0.385579371 0.297777778
## beta[98] -2.642365e-01 -0.83566182 0.247318327 0.201111111
## beta[99] -3.475842e-01 -0.89264096 0.177701067 0.118888889
## beta[100] -4.208126e-01 -1.04027357 0.144216541 0.076666667
## beta[101] -4.899116e-01 -1.14609123 0.067911356 0.050000000
## beta[102] -5.653700e-01 -1.25525846 0.041526380 0.031111111
## beta[103] -6.329116e-01 -1.36135924 -0.005417412 0.024444444
## beta[104] -6.894248e-01 -1.44018748 -0.051513452 0.016666667
## beta[105] -7.259733e-01 -1.49806783 -0.059065080 0.011111111
## beta[106] -7.467424e-01 -1.44497269 -0.044509144 0.013333333
## beta[107] -7.537286e-01 -1.47905222 -0.037769138 0.016666667
## beta[108] -7.792557e-01 -1.55581617 -0.018457969 0.020000000
## beta[109] -8.001735e-01 -1.55813086 -0.054858819 0.014444444
## beta[110] -8.113561e-01 -1.62778758 -0.029250561 0.016666667
## beta[111] -8.337856e-01 -1.64641481 -0.049182394 0.016666667
## beta[112] -8.442536e-01 -1.81183085 -0.006933733 0.024444444
hist(ej6c.bugs$sims.list$beta)
Ejercicio 6d
cat Ej6d.txt
## model
## {
## #Likelihood
## #Space eq.
## for (i in 1:n) {
## y[i] ~ dpois(mu[i])
## }
## #State eq.
## for (i in 2:n) {
## mu[i] ~ dgamma(tau.b,b[i])
## b[i] <- tau.b/mu[i-1]
## }
## #Priors
## mu[1] ~ dgamma(0.001,0.001)
## #tau.b ~ dgamma(0.001,0.001)
## tau.b ~ dgamma(200,1)
##
## #Prediction 1
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
##
## }
#-Defining inits-
# inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}
#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
parameters<-c("tau.b","yf1","mu")
ej6d.bugs<-bugs(data,inits,parameters,model.file="Ej6d.txt",
n.iter=5000,n.chains=1,n.burnin=500)
ej6d.jags<-jags(data,inits,parameters,model.file="Ej6d.txt",
n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "x" in data
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 112
## Unobserved stochastic nodes: 225
## Total graph size: 453
##
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
# print("Bugs")
# exploracionMCMC(ej6d.bugs)
# print("Jags")
# exploracionMCMC(ej6d.jags)
Betas
# betasMCMC(ej6d.bugs)
# betasMCMC(ej6d.jags)
DIC
dicMCMC(ej6d.bugs)
## [1] 343.8262
dicMCMC(ej6d.jags)
## [1] 342.861
# graficarY1(ej6d.bugs)
# graficarY1(ej6d.jags)
DICs
dicMCMC(ej6a.bugs)
## [1] 347.5122
dicMCMC(ej6a.jags)
## [1] 399.8575
dicMCMC(ej6b.bugs)
## [1] 345.9282
dicMCMC(ej6b.jags)
## [1] 346.0538
dicMCMC(ej6c.bugs)
## [1] 338
dicMCMC(ej6c.jags)
## [1] 346.4019
dicMCMC(ej6d.bugs)
## [1] 343.8262
dicMCMC(ej6d.jags)
## [1] 342.861
